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The  research  described  in  this  report  was  the  continuation  of  an 
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Group  headed  by  David  K.  Bowser.  This  work  was  performed  in  the  Depart¬ 
ment  of  Electrical  Engineering,  University  of  Kentucky,  under  Air  Force 
Grant  AFOSR-82-0207.  This  report  covert;  work  performed  from  June  1, 
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May  1,  1982  to  May  15,  1982. 
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digital  controller  was  written  by  Dr.  Kyung-Sam  Choi,  Associate  Professor 
of  Electrical  Engineering,  Hong-ik  University,  Republic  of  Korea,  when 
he  was  visiting  the  University  of  Kentucky  in  the  academic  year  1980-1981. 
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SECTION  I 


\ 

\  INTRODUCTION 

j 

The  digitalization  of  flight  control  systems  has  been  of  increasing 
interest  to  the  Air  Force.  One  of  the  problems  confronting  the  designer 
is  the  real-time  implementation  of  advanced  control  algorithms  within 
the  computational  capability  of  the  on-board  computer.-"  Although  digital 
^computer'  technology  continues  to  advance  significantly,  new  and  expanded 
software  requirements  for  such  functions  as  navigation,  display  and 
control  manage  to  keep  pace  with  Improvements  in  computational  capa¬ 
bility.  As  a  rule,  only  a  small  fraction  of  the  CPU  frame  time  is 
allocated  for  control  law  computation.  Hence,  from  the  standpoint  of 
implementation,  the  sampling  rate  should  be  sufficiently  low  in  order 
to  allow  time  for  computation  and  for  the  computer  to  be  time -shared. 

But  lower  limits  of  the  sampling  rate  are  determined  by  factors  such  as 
excessive  ratcheting  and  jittering  in  the  time  response,  errors  due  to 
measurement  noise  and  sensitivity  to  plant  parameter  uncertainty  and 
disturbances.  The  con’trol  system  designer  is  then  faced  with  the  task 
of  optimizing  the  performance  of  the  digital  control  system  at  a  given 
rate  of  sampling. 

""  *In  converting  a  continuous-data  (analog)  controller  into  a  digital 
controller,  ad  hoc  approaches  such  as  bilinear  transform  and  prewarped 
Tustin  transform  techniques  have  typically  been  used.  These  methods 
have  the  advantage  of  being  straightforward  and  easy  to  use,  and  they 
are  intuitively  appealing.  But  the  performance  of  a  system  digitalized 
by  these  approaches  resembles  the  performance  of  the  baseline  (continu¬ 
ous)  system  only  when  the  sampling  frequency  is  relatively  high, 
because  the  dynamics  of  the  plant  and  the  feedback  structure  of  the 
system  are  not  taken  into  consideration.  _ _ 

Four  years  ago  Rattan  and  the  author  of  this  report  presented  a 
method  [1]  using  a  complex-curve  fitting  technique  to  synthesize  the 
digital  controller  so  that  the  frequency  response  of  the  digitalized 
system  matches  that  of  the  original  continuous  model  with  a  least-square 
fit.  This  method  was  recently  applied  to  a  mathematical  model  of  the 
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longitudinal  flight  control  of  the  YF-16  fighter  aircraft  at  an  alti¬ 
tude  of  30,000  ft.  and  Mach  0.6.  Results  [2]  better  than  that  of  Tustin 
transform  approach  have  been  obtained,  especially  for  lower  sampling 
frequencies.  However,  this  method  does  not  take  the  time-domain  perfor¬ 
mances  into  consideration;  and  only  the  magnitude  plots  of  the  frequency 
responses  are  matched,  without  regard  to  the  phase  plot.  Moreover,  the 
comparison  between  the  frequency  response  of  a  continuous  system  and  the 
frequency  response  of  a  digital  control  system  becomes  meaningless  (see 
Appendix  D)  as  the  signal  frequency  approaches  the  folding  frequency 
(one  half  of  the  sampling  frequency) .  To  compensate  for  these  short¬ 
comings,  the  state-variable  design  techniques  in  the  time-domain  should 
be  developed. 

In  a  previous  study,  the  author  and  his  associate  showed  that  [3] 
by  using  z-transformation  technique,  a  digital  controller  can  always  be 
synthesized  so  that  the  output  of  the  digitalized  system  matches  the 
output  of  the  continuous  model  at  all  sampling  instants,  under  the 
stimulation  of  the  same  input. 

In  this  study,  an  attempt  is  made  to  match  the  continuous  state 
trajectory  of  the  digital  control  system  with  that  of  its  analog 
(continuous-data)  model.  Matching. the  state  trajectories  instead  of 
the  output  responses  assures  that  the  performances  of  the  internal 
variables  of  the  plant,  as  well  as  the  output  variable,  are  preserved 
in  the  digitalization.  It  should  also  be  emphasized  that  the  matching 
is  specified  over  the  entire  continuous  time  axis,  not  just  at  discrete 
sampling  instants,  and  is  quantified  by  a  minimum  integral  squared 
error.  The  choice  of  this  performance  criterion  is  motivated  by  the 
fact  that  if  the  state  trajectories  of  two  linear  dynamical  systems 
match,  then  frequency  responses  of  the  two  systems  will  also  match,  as 
seen  by  Laplace-transforming  the  state  equations. 

The  mathematical  tool  used  in  this  research  is  an  extended  maximum 
principle  of  the  Pontryagln  type,  which  enables  one  to  synthesize  a 
"staircase"  type  of  optimal  control  signals,  such  as  the  output  signal 
of  a  zero-order  hold  associated  with  a  digital  controller.  The  extended 
maximum  principle  was  initiated  by  Chang  [4]  and  further  developed  by 
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the  author  and  his  co-workers  [5-9],  Though  not  widely  used  in  the 
technical  literature,  it  is  a  very  effective  tool  for  obtaining  the 
discrete  time  optimal  control  policy  for  an  integral  performance  index, 
with  a  minimum  amount  of  mathematical  manipulations.  It  is  shown  in 
the  sequel  (Section  VI)  that  the  same  result  can  be  derived  from  the 
discrete  maximum  principle  and  an  equivalent  discrete-time  performance 
index  but  through  a  far  more  cumbersome  procedure. 

This  research  is  a  continuation  of  the  effort  initiated  in  the  1981 
Summer  Faculty  Research  Program  sponsored  by  the  Air  Force  Office  of 
Scientific  Research  through  the  Southeastern  Center  for  Electrical 
Engineering  Education,  under  contract  F49620-79-C-0038.  In  order  for 
this  report  to  be  self  contained,  the  mathematical  formulation  of  the 
problem  will  be  repeated  (Section  II).  The  optimal  strategy  and  the 
digital  controller  are  rederived  for  the  general  case  (Sections  III, 

IV).  Computational  verification  and  performance  evaluations  and  com¬ 
parisons  are  presented  in  Section  V,  where  computational  difficulties 
that  are  special  to  this  method  are  discussed.  Conclusions  and  recom¬ 
mendations  are  given  in  Section  VII. 
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SECTION  II 


FORMULATION  OF  THE  DESIGN  PROBLEM 


Consider  a  continuous  output  tracking  system  (Fig.  1)  that  has 
satisfactory  (or  ideal)  performances.  The  state  and  output  equations 


of  the  plant  are  given  by 

x  (t)  =  Ax  (t)  +  bu  (t)  (1) 

-a  -a  —  m 

y  (t)  =  cx  (t)  +  du  (t)  (2) 

m  —a  m 

The  state  and  output  equations  of  the  controller  are  given  by 

xc(t>  =  A^t)  +  £cem(t)  (3) 

um(t)  =  ZcZc(t)  +  dcem(t)  (4) 


where  x  (t)  and  x  (t)  are  n  and  n  dimensional  vectors,  respectively, 
-a  -c  c 

and  u  (t),  e  (t),  y  (t)  and  r(t)  are  scalar  functions.  The  dimensions 
m  m  m 

of  the  coefficient  matrices  are  commensurate  with  the  vectors  with 
which  they  associate.  The  design  objective  is  to  replace  the  control¬ 
ler  G  (s)  by  a  digital  controller  D(z)  such  that  the  state  trajectory 
c 

of  the  digitalized  system  matches  that  of  the  continuous  model  as 
closely  as  possible.  The  digital  control  system  is  represented  by 
Fig.  2,  where  G(s)  is  the  same  plant  as  in  the  continuous  model,  and 
D(z)  is  to  be  synthesized  in  such  a  way  that  when  r(t)  is  a  unit-step 
function,  the  performance  index 

00 

J  -  4/  {[X(t>  -  x  (t)]'Q[x(t)  -  x(t)] 

i  Q  -  -a  -  -a 

+  B[u(t)  -  u(t)]2}  dt  (5) 

m 

attains  its  minimum,  where 


4 


r  (t) 


*c(t) 


X  (t) 
—o 


e  (t ) 
m 


Tr  1 


6c(s) 


u  (t) 
m 


G  (s) 


y  (t) 

m 


continuous 

controller 


plant 
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Fig.  2  The  Digitalized  System 
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Note  that  the  performance  index  is  an  integral,  not  a  discrete  sum. 
Therefore  an  attempt  is  made  to  match  the  trajectories  over  the  con¬ 
tinuous  time  axis,  not  just  at  the  sampling  instants.  The  state  and 
output  equations  of  the  plant  in  the  digital  control  system  are 

x(t)  -  Ax(t)  +  bu(kT)  (7) 

y(t)  »  cx(t)  +  du(kT)  (8) 

for  kT  <  t  <  (k  +  1)T,  on  account  of  the  zero-order  hold  used  in  the 
digital  control  system  (Fig.  2). 
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SECTION  III 


THE  OPTIMAL  CONTROL  SEQUENCE 


(A)  The  Extended  Maximum  Principle 

An  extended  version  of  the  maximum  principle  of  Pontryagin  will 
be  used  to  find  the  optimal  control  sequence  u(kT),  k  -  0,1,2,*«», 
which  minimizes  the  performance  index  (5).  The  error  sequence  e(kT) 
can  be  expressed  in  terms  of  u(kT),  and  the  digital  controller  D(z)  can 
be  determined  by 

D<2>  •  w 

where  U(z)  is  the  z-transform  of  u(kT)  and  E(z)  is  the  z-transform  of 
e(kT). 

The  original  Pontryagin’ s  Maximum  Principle  cannot  be  applied  to 
this  problem  in  which  the  control  function  must  be  sampled  and  held, 
i.e. , 

u(t)  -  u(kT)  ,  kT  <  t  <  (k  +  1)T  (10) 

The  derivation  of  Pontryagin' s  maximum  principle  is  based  on  the 
assumption  that  piece-by-piece  patchwork  of  admissible  control  func¬ 
tions  is  again  admissible.  This  condition  is  not  satisfied  by  the 
output  of  a  zero-order  hold  following  a  constant  rate  sampler,  if  the 
point  of  patching  is  selected  at  an  instant  other  than  the  sampling 
instants.  However,  the  extended  maximum  principle  may  be  applied  to 
the  case  where  the  control  inputs  are  outputs  of  zero-order  holds 
[5-9J.  The  following  derivation  follows  from  Reference  [7].  The 
Hamiltonian  function  H  is  formed  from  the  adjoint  vector  j>(t),  the 
state  derivative  x(t)  and  the  integrand  of  the  performance  index  as 

H(x(t),£(t),u(t))  -  £’(t)x(t) 

-  -|{(x(t)  -  xa(t)l’Q(x(t)  -  xa(t)J  +  6(u(t)  -  uffl(t)]2}  (ID 
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where  the  prime  denotes  transpose.  Substituting  (7)  into  (11)  gives 


H(x(t) ,p(t) ,u(t))  -  p'(t)Ax(t)  +  £*bu(t) 

-  4(tx(t)  -  x  (t)]’Q[x(t)  -  x  (t)]  +  6[u(t)  -  xi At)}2}  (12) 

l  —  —a  —  —a  m 

The  adjoint  vector  g(t)  is  given  by 

2(b)  B  H(x(t)  ,£(t)  ,u(t))  (13) 

Substituting  (12)  into  (13)  gives  the  adjoint  system 

p(t)  -  -A'£(t)  +  Q[x(t)  -  xa(t)]  (14) 

The  necessary  condition  for  the  control  and  trajectory  to  be  optimum  is 
that  the  first  variation  of  the  integral  of  the  Hamiltonian  is  zero, 
i.  e . , 

l  “ 

lim  —  /  [H(x,g,u  +  <5  u)  -  H(x,g,u)]  dt  *  0  (15) 

e-K)  E  0  E 

where  the  variable  t  in  the  Hamiltonian  has  been  suppressed,  and 

6£u(t)  *  eAu(kT)  ,  kT<t<(k  +  l)T  (16) 

where  Au(kT)  is  an  arbitrary  finite  variation  of  u(kT).  Substituting 
(12)  and  (16)  into  (15)  gives 

«  (k+l)T 

I  /  l£'b  -  6(u(kT)  -  u  (t))]Au(kT)  dt  -  0  (17) 

k-0  kT 

Since  Au(kT)  is  an  arbitrary  constant  which  may  be  set  independently  for 
each  k,  (17)  implies 
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(k+l)T 

/  [g'b  -  6(u(kT)  -  u  (t))]  dt 

kT 


0 


(18) 


From  (18)  Che  optimal  control  u(kT)  is  found  to  be 


,  (k+l)T 

u(kT)  -  -£=  I  [p'b  +  8u  (t) ]  dt  (19) 

eT  kT 


The  difficulty  in  the  application  of  the  maximum  principle  is  the  two- 
point  boundary  value  problem.  In  this  formulation,  the  two-point 
boundary  value  problem  is  as  follows:  For  any  i,  if  x^(0)  is  given, 
then  p^(0)  is  unspecified  (remains  to  be  determined  from  (7),  (19)  and 
(14)).  If  xi(0)  is  unspecified  (remains  to  be  determined  by  the 
maximum  principle),  then  pi(0)  ■  0.  The  same  rule  applies  to  x^(«)  and 
P1(°’). 

(B)  Determination  of  the  Optimal  Control  Sequence 


In  order  to  compute  the  optimal  control  sequence  u(kT)  from  (19) , 
the  solution  p(t)  of  (14)  and  an  explicit  expression  of  um(t)  must  be 
found  first. 


Let  the  augmented  state  vector  of  the  model  be 


x(t) 


*a(t) 


Sc(t) 


(20) 


It  follows  from  Appendix  A  that 


i  (t)  -  A  x  (t)  +  b  r(t)  (21) 

•in  ni“iu  “in 

u  (t)  ■  c  x  (t)  +  d  r(t)  (22) 

m  nr— in  to 


where 


9 


and 


d 


m 


d 

c 

1+dd 

c 


Define  Q  to  be  the  n*(n  +  n  )  matrix  obtained  by  augmenting  n 
me 

of  zeros  to  Q,  i.e., 

Qg  *  IQ  °J 

Then  by  definition 

Qx  (t)  -  Q  x (t) 

—a  m— m 


Assume  that  r(t)  is  a  step  function,  i.e.. 


r(t) 


/  a  t  >0 

(  0  t  <0 


It  is  shown  in  Appendix  B  that,  for  kT  <  t  $  (k  +  1)T, 


(23) 


(24) 


(25) 


(26) 


columns 

c 


(27) 


(28) 


(29) 
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x(t)  -  $(t-kT)x(kT)  +  $  (t-kT)bu(kT) 


x  (t)  *  <t  (t-kT)x  (kT)  +  $  (t-kT)b  a  (31) 

—in  m  —ni  in  — m 


where 


¥'(t-kT)p(kT)  +  F(t-kT)x(kT)  +  F  (t-kT)bu(kT) 


-  F  (t-kT)x  (kT)  -  F  (t-kT)b  a 
in  “in  m  “in 


$(t)  =  e 


*S(t) 


*(t) 

m 


I  $(t)  dT 
0 


4>S(t)  -  /  $  (t)  dr 


y(t)  -  e 


FS(t) 


F  (t) 
n 


F®(t) 

in 


/  ’f(t-x)Q#(x)  dx 

0 


/  F(t)  dT 
0 


/  'F(t-T)Q  4>  (t)  dx 
0  m  m 


/  F  (t)  dx 
0  m 


*(t)  *Q«(t) 


▼(t)  *Q  *  <*) 

m  m 


Substituting  (31)  and  (32)  into  (19)  and  solving  for  u(kT)  yields 


u(kT)  -  - - - { b,'fS(T)p(kT)  +  b’FS(T)x(kT) 

BT  -  b’FSS(T)b  " 


+  [Be  *S(T)  -  b‘ FS(T)  ]x  (kT) 

in  m  *•  uj  —in 


+  [Be  $SS(T)b  +  BTd  -  b'FSS(T)b  ]al  (42) 

mm  -m  m  —  m  — m  » 


where 


4,S(T) 

T 

«  / 

f(t)  dt 

(43) 

0 

FSS(T) 

T 

“  / 

Fs(t)  dt 

(44) 

0 

FaS(T) 

m 

T 

Fa(t)  dt 
m 

»  / 
n 

(45) 

and 


*SS(T)  -  /T  *S(t)  dt 

(46) 

m  ;0  m 

Equation  (42)  may  be  used  to  compute  u(kT)  iteratively  after  the  two- 
point  boundary  problem  is  solved.  In  this  research,  (9)  is  used  to 
synthesize  the  digital  controller  transfer  function  D(z).  The  two- 
point  boundary  problem  still  arises  and  will  be  treated  in  the  next 
section. 
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SECTION  IV 


THE  DIGITAL  CONTROLLER 


(A)  z-Transform  of  the  Optimal  Control  Sequence 

In  order  to  determine  U(z)  for  use  in  Eq.  (9),  P(z),  X(z)  and 
Xm(z)  must  be  determined  first.  Setting  t  *  (k  +  1)T  in  (30)  -  (32) 
gives 


x[(k+l)T]  -  4>(T)x(kT)  +  4>  (T)bu(kT) 


x  t(k+l)T]  -  «  (T)x  (kT)  +  4>  (T)b  a 

-m  m  -m  m  —in 


£[(k+l)T]  -  4<(T)£(kT)  +  F(T)x(kT)  -  Fffl(T)x  (kT) 


+  FS(T)bu(kT)  -  FS(T)b  a 
—  m  — m 


Taking  z- transformation  of  (47)  -  (49)  gives,  respectively, 
X(z)  -  $(z)[zx(0)  +  *S(T)bU(z)] 


X  (z)  *  * 


P(z) 


'  (z)  fzx  (0)  +  <(T)b 

m  l  m  — m  z— l  I 

?(z)  jj£(0)  +  F(T)X(z)  -  Fm(T)Xin(z) 

az  **1 
z-lj 


+  FS(T)bU(z)  -  FS(T)b 

-  m  — m 


where 


(47) 

(48) 


(49) 


(50) 

(51) 


(52) 


$(z)  * 

[ zl  -  *(T)]"1 

(53) 

-e-  > 

B 

/-s 

N 

V 

[zl  -  ♦(T)]“1 

m 

(54) 

-€> 

N 

1 

[zl  -  KT)]-1 

(55) 
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Substituting  (50)  and  (51)  into  Eq.  (52)  gives 


P(z)  -  'P(z){zp(0)  +  F(T)$(z)zx(0)  -  F  (T)4>  (z)zx  (0) 

—  •  **  m  m  -m 

-  lFm (D* -(z)*®(T)  +  F®(T)]b 

in  rn  m  m  "n  zx 

+  lF(T)*(z)4>S(T)  +  FS(T)  ]bU(z) }  (56) 

Substituting  (50),  (51)  and  (56)  into  the  z-transform  of  (42),  and 
solving  for  U(z)  yields 

U(z)  "  gf^pH(I)b  {b,K(z)zx(°)  +  brK0(z)zE(0) 

+  [6c  4>8(T)$  (z)  -  b’K  (z)]zx  (0) 
mm  m  ~  m  —m 

+  [gc($S(T)$(z)*®(T)  +  ^S(T))bm  -  b’H  (z)b  +  BTd  ]  (57) 

mm  m  m  m  — m  -  m  — m  m  z— i 


where 

K  (z)  -  [¥S(T)?(z)F  (T)  +  F®(T)]J  (Z)  (58) 

m  m  m  m 

K(z)  -  [1*S(T)?(z)F(T)  +  FS(T)]S(z)  (59) 

Kq(z)  -  4,S(T)?(z)  (60) 

H  (z)  "  f!S(T)  +  Kn(z)F®(T)  +  K  (z)*"(T)  (61) 

m  m  u  m  mm 

H(z)  -  FSS(T)  +  K0(z)FS(T)  +  K(z)*®(T)  (62) 

(B)  z-Transform  of  the  Error  Sequence 


From  the  block  diagram  of  Fig.  2  and  Eqs.  (8)  and  (30)  we  may 

write 
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e(t)  -  r(t)  -  y(t)  -  ct  -  y(t)  (63) 

y(t)  *  c$(t-kT)x(kT)  +  [c$S(t-kT)b  +  d]u(kT)  (64) 

for  kT  <  t  <  (k  +  1)T.  Therefore 

e[(k+l)T]  -  a  -  c$(T)x(kT)  -  [c«S(T)b  +  d]u(kT)  (65) 

Taking  z-transform  of  the  above  equation  gives 

E(z)  -  e(0+)  +  -At  -  -  c*(T)X(z)  -  -  [c*S(T)b  +  d]U(z)  (66) 

Z“ 1  z  z  ~ 


where  e(0+)  is  found  from  (63)  and  (8)  as 

e(0+)  -  a  -  cx(0+)  -  du(0+)  (67) 

Substituting  (67)  and  (50)  into  (66)  gives 

E(z)  -  -  du(0+)  -  c[I  +  $(T)J(z)]x(0+) 

Z“1  — 

-  -tc(I  +  $(T)*(z))4>S(T)b  +  d]U(z)  (68) 

z  - 

From  the  definition  of  $(z)  given  in  (53),  it  follows  that 

I  +  4>(T)$(z)  ■  z$(z)  (69) 

Moreover,  for  most  control  systems  the  control  signal  does  not  supply 
power  to  the  output.  Hence  d  *  0  and  E(z)  of  (68)  may  be  written  as 

E(z)  -  -  cz$(z)x(0+)  -  c$(z)*S(T)bU(z)  (70) 

z— i  —  — 


By  virtue  of  the  relationship  given  in  (9),  the  z-transfer  function  of 
the  digital  controller  may  now  be  written,  provided  that  the  initial 
conditions  x(0),  xffl(0)  and  g(0)  are  known. 
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(c)  The  z-Transfer  Function  of  the  Digital  Controller 


I 


In  the  model-following  design,  the  initial  state  x  (0)  of  the  model 
may  be  taken  to  be  zero.  But  if  x(0)  is  given,  the  p(0)  cannot  be 
arbitrarily  specified.  The  maximum  principle  requires  that  p(0)  be  such 
that  Eqs.  (7),  (14)  and  (19)  yield  an  optimal  control  sequence  that 
drives  the  system  from  the  given  x(0)  to  a  prespecified  final  state 
x(«) ,  or  results  in  £(«)  ■  0  if  the  final  state  is  unspecified  (free). 
Numerical  determination  of  g(0)  for  given  x(«)  or  g(«)  is  virtually 
impossible  for  systems  requiring  a  large  number  of  sampling  periods  to 
run,  because  either  matrix  A  or  -A*  has  eigenvalues  in  the  unstable 
region.  On  the  other  hand,  if  x(0)  is  unspecified,  g(0)  must  be  set 
equal  to  zero.  But  then  x(0)  is  determined  by  Eqs.  (7),  (14)  and  (19) 
and  the  final  condition,  which  is  either  at  a  given  value  of  x(«°)  or  at 
£(a>)  =  o.  Again,  numerical  determination  of  x(0)  is  impractical  for 
systems  requiring  a  large  number  of  sampling  periods  to  run. 

For  the  system  under  consideration,  the  state  trajectory  of  the 

continuous  model  is  assumed  to  start  from  x  (0)  ■  0.  Let  x(0),  the 

-m  — 

initial  state  of  the  plant  of  the  digital  control  system,  be  unspeci¬ 
fied.  Then  g(0)  ■  0.  However,  since  the  minimization  of  the  perfor¬ 
mance  Index  given  in  (5)  means  the  continuous  matching  of  x(t)  with 
x  (t)  over  an  infinitely  long  period  of  time,  it  is  reasonable  to  con- 

” O 

jecture  that  x(t)  starts  at  the  same  point  as  x  (t),  or  very  close  to 

•  “fl 

it,  provided  that  the  sampling  frequency  is  considerably  higher  than  the 
natural  frequency  of  the  control  system.  Hence  the  initial  condition 
may  be  chosen  as 

x  (0)  -  0  (71) 


x(0)  -  0 


E(0)  -  0 


Substituting  these  conditions  into  Eqs.  (57)  and  (70),  and  using  the 
resulting  expressions  in  the  right-hand  side  of  (9),  we  obtain 


16 


I 


D(z) 


E(z) 

U(z) 


-c$(z)$  (T)b 


ST  -  b’H(z)b 


8cm[$®(T)S(z)$®(T)  +  $!S(T)]b 
in  m  in  in  m  —in 


b’H  (z)b  +  STd 
—  in  — m  in 


(74) 


It  should  also  be  recognized  that 


c$(z)4>S(T)b  +  d  =  z-transform  of  |  — % - G(s)  |  (75) 


G<s)] 


For  S  “  0,  the  digital  controller  is  obtained  from  (74)  as 


D(z) 


b'H(z)b*  I  " 

iiw  - '♦<*>♦  (T)i 

-  m  -m 


(76) 


Note  that  Eq.  (19)  cannot  be  used  to  find  u(kT)  for  8  =  0.  But  Eq.  (18) 
becomes 

(k+l)T 

/  p’bdt  =  0  (77) 

kT 

Since  p  depends  on  u(kT) ,  (77)  can  be  used  to  derive  the  optimal  control 
sequence.  This  approach  was  taken  in  the  previous  project  [10]  and  the 
digital  controller  of  (76)  was  obtained. 


(D) 


The  Closed-Loop  Transfer  Function 
Let  G^G(z)  denote  the  z-transform  of 


G(s), 


i.e.. 


Let 


GhG(z) 


c4>(z)$S(T)b 


(78) 


G  (z) 
m 


b’H  (z)b 
-  m  -m 

b*H(z)b 


(79) 
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The  z-transfer  function  of  the  closed-loop  system  of  Fig.  2  is  given  by 


Y(z)  _  P(z)GhG(z) 

R(z)  1  +  D(z)gTg(z) 
n 


(80) 


If  the  digital  controller  of  (76)  is  used,  D(z)  may  be  written  as 


D(z) 


G(z) 

_  m _ 

1  -  G  (z)gTg(z) 
m  n 


Substituting  (81)  into  (80)  yields 


(81) 


Y(z) 

R(z) 


G  (z)G.G(z) 
m  n 


b'H  (z)b 


—  m 


b'H(z)b 


GhG(z) 


(82) 


Thus  the  optimal  model-following  design  of  digital  controller  can  be 
achieved  by  the  Guillemin-Truxal  approach  [11]  which  selects  a  digital 
controller  transfer  function  algebraically  so  that  the  closed-loop  z- 
transfer  function  is  the  product  of  the  z-transfer  function  of  the  plant 
with  a  zero-order  hold  and  a  correction  factor  Gffi(z)  given  by  (79). 

It  is  readily  observed  from  (58)  -  (62)  that 

_A’t  at 

(1)  Poles  of  H  (z)  are  eigenvalues  of  e  and  e  111  , 

m 

-A'T  AT 

(2)  Poles  of  H(z)  are  eigenvalues  of  e  and  e 

-  AT 

Since  the  poles  of  G^G(z)  are  eigenvalues  of  e  ,  what  the  multiplication 

factor  G  (z)  in  (82)  actually  does  is,  among  other  things,  to  replace  the 
m _ 

poles  of  G^G(z)  by  z-plane  poles  of  the  continuous  model  with  samplers 
connected  to  in  input  and  the  output.  Thus  the  closed-loop  digitalized 
system  will  have  poles  of  the  sampled  continuous  model,  which  are  eigen¬ 
values  of 

The  stability  conditions  of  the  digitalized  system  is  given  in  [12]. 
The  system  is  easily  stabilized  if  all  diagonal  elements  of  Q  are  positive. 
It  is  also  observed  from  (82)  that  if  the  closed-loop  system  is  stable, 
then  any  zeros  of  b'H(z)b  outside  the  unit  circle  must  be  cancelled  by 


18 


some  zeros  of  b'H  (z)b  .  Theoretically,  zeros  of  b'H(z)b  outside  the 
—  m  — m  _ —  — 

unit  circle  may  be  cancelled  by  some  zeros  of  G,G(z )  but  such  a  system  is 

h 

not  realizable. 


(E)  Controller  In  the  Feedback.  Path 

When  the  controller  is  in  the  feedback  path,  as  in  Fig.  3  and  Fig.  4, 
Formulas  (74)  and  (76)  no  longer  apply  and  should  be  rederived.  In  this 
case,  the  system  equations  (1)  -  (4)  and  (7)  -  (8)  are  replaced  by 


and 


“  A  x  (t)  +  b  e  (t) 

— 3.  —  m 

(83) 

=  c  x  (t)  +  d  e  (t) 

— a  in 

(84) 

jr 

■  A  x  (t)  +  b  y  (t) 
c  — c  — c  m 

(85) 

V  (t) 
m 

=  c  x  (t)  +  d  y  (t) 
c  — c  c  m 

(86) 

£(t)  1 

■  A  x(t)  +  b  e(kT) 

(87) 

y(t)  < 

■  c  x(t)  +  d  e(kT) 

(88) 

respectively.  The  derivation  follows  a  similar  pattern  as  that  of  (74) 
and  (76),  and  will  not  be  repeated  here.  The  result  for  the  case  of 
8  *  0  is  given  as  follows: 


D1(z) 


Vizi 

Y(z) 


b'H  (z)b 
l  ~  m  ~m 

b’H(z)b 


z  -  1 


b’F®S(T)b 
—  m _  a 

b'FSS(T)b 


_  .  b’H  (z)b 

d  +  tcw?(T)b.i] 


Z  -  1 


(89) 


Again,  for  d  ■  0, 
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controller 


(90) 


D^z) 


b' [H(z)b  -  H  (z)b  ] 
hi  m 

tc  *  (z)*s(t)b] [b'H  (z)b  ] 
—  in  tn 


where  H(z),  H  (z),  $(z),  and  *S(T)  are  again  given  by  (62),  (61),  (53), 
m 

and  (34),  respectively,  but  A  and  b  are  now  given  by: 

n  ~"in 


A 


m 


b 

— m 


bd  c 

A  -  - 

1  +  dd 

c 

b  c 
c 

1  +  dd 

c 

b 

1  +  dd 

c 

b  d 
c 

1  +  dd 

c 


-be 
—  c 

1  +  dd 

c 

b  dc 
~c  c 

“  1  +  dd 


(91) 


(92) 


which  are  different  from  their  previous  definitions  of  (23)  and  (24). 

It  is  interesting  to  see  that  for  the  control  system  of  Fig.  4  with 
a  digital  controller  given  by  (90) ,  the  closed- loop  transfer  function  is 
again  given  by  (82). 

Careful  examination  of  (90),  in  view  of  (58)  -  (62),  shows  that  the 
digital  controller  D^(z)  is  not  realizable,  namely,  the  numerator  poly¬ 
nomial  of  D^(z)  is  one  degree  higher  than  that  of  the  denominator.  To 
remedy  for  this,  the  digitalized  system  can  be  implemented  as  in  Fig.  5, 
which  can  be  derived  as  follows.  First,  V(z)  of  Fig.  4  is  found  to  be 


V(z) 


b'H  (z)b 

1  —  m  — m  otz 

b'H(z)b  z  -  1 


(93) 


In  view  of  (58)  -  (62),  the  initial  value  v(0)  is  found  to  be 


v(0)  »  V(z) 

Z 


L  b’FSS(T)bJ 


(94) 
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D,  (z) 


Fig.5  Realization  of  The  Feedback  Digital  Controller  D(  (z) 
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Thus,  when  r(t)  is  a  step  function  of  magnitude  a 


D  2(z) 


V(z)  - 

z  —  X 


Y(z) 


b'FSS(T)b 
—  m  — m 

b'F(T)b 


b’H(z)b 

b '  H(  z  )  b 
—  m 


-  1 


c<t(z)$8(T)b 


(95) 


The  digital  controller  D2U )  of  (95)  is  realizable,  on  account  of  (53)  - 
(55)  and  (58)  -  (62). 


SECTION  V 


COMPUTATIONAL  VERIFICATION 


Attempts  have  been  made  to  verify  the  end  results  through  the 
numerical  computations  of  the  transfer  functions  of  the  digital  con¬ 
trollers  and  the  evaluations  of  the  performances  for  three  examples. 

For  the  two  lower  order  systems,  excellent  time  responses  of  the  overall 
closed-loop  digital  control  systems  and  satisfactory  frequency  responses 
were  obtained.  The  other  example  was  the  inner  loop  of  the  longitudinal 
control  of  YF-16  fighter  aircraft  at  30,000  ft  altitude  at  Mach  0.6. 
Limited  by  the  capability  of  the  available  computer  programs,  the  re¬ 
sult  for  the  last  example  was  unsatisfactory,  hence  not  presented.  The 
computational  difficulties  and  recommendations  will  be  discussed  later. 

In  comparing  the  frequency  responses  for  the  following  examples, 
the  expression 


Me1"1) 


D(z)G,G(z) 

1  +  D(z)gTg(z) 


Z-*£ 


ju>T 


(96) 


is  used  for  the  frequency  response  of  a  digital  control  system  represented 
by  Fig.  2.  It  should  be  noted  that,  in  light  of  the  concept  developed 
in  Appendix  D,  there  does  not  exist  a  closed-form  expression  of  the  fre¬ 
quency  response  of  a  digital  control  system.  It  is  neither  represented 
by  (96),  nor  is  it  represented  by 


IM 

R(jm) 


D(eJ«T)  -  G(j») 

1  +  D(cJ“T)GhG(ej“T) 


(97) 


Example  1.  The  end  result  given  in  (74)  or  (76)  is  computationally  com¬ 
plicated.  It  will  be  demonstrated  through  simple  examples  first.  Refer¬ 
ring  to  Fig.  1  and  Fig.  2,  let  the  plant  and  the  controller  be 

GO  -  y-J-J  (98) 
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and 


(99) 


G  (s) 
c 


8  +  2 
s 


The  coefficients  in  the  state  variable  model  are 

c  -  [1]  d  *  0  (100) 

c  »  [21  d  -  1  (101) 

c  c 

We  choose  Q  =  1,  ■  [1  0]  and  8  *  0.  The  sampling  period  is  chosen 

to  be  relatively  large  because  the  system  is  a  slow  one.  Otherwise, 
conversion  to  digital  control  would  not  have  been  a  problem.  We  shall 
consider  the  cases  of  T  *  0.5  sec.,  T  ■  1  sec.  and  T  =  2  sec. 

Case  1.  T  *  0.5  sec. 

Implementation  of  (76)  on  the  computer  program  gives 

n .  .  1. 1017 (z  -  0.3682)(z+  0.3078)  (z  -  0.6065)  (2  +  4.1231)  , 

U{Z)  ”  (z  -  l)(z  +  0. 2769)  (z  -  0. 5399)  (z  +  3.8149) 

Since  the  closed-loop  system  is  stable  (integral-squared  error  with  re¬ 
spect  to  a  stable  model'  is  minimized)  ,  the  pole  and  zero  of  (102)  out¬ 
side  the  unit  circle  should  not  be  implemented.  Since  the  performance 
of  the  digital  control  system  should  match  its  continuous  model  only  at 
low  frequencies,  the  controller  transfer  function  of  (102)  may  be  approx¬ 
imated  by  setting  z  ■  1  in  the  factors  (z  +  4.1231)  and  (z  +  3.8149). 

Thus 


A  -  [-1]  b  -  [1] 

Ac  -  [0]  ^  -  [1] 


,  .  1. 1722 (z  -  0, 3682) (z  +  0.3078)(z  -  0.6065)  , 

DW  "  (z  -  1)  (z  +  0. 2769)  (z  -  0.5399)  k 

For  the  purpose  of  comparison,  the  controller  G  (s)  is  discretized  by 

c  2  z  -  1 

prewarped  Tustin  transformation  (substituting  s  by  ^  ^  and  2  by 

T 

2  x  2  respectively,  in  Gc(s)): 
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(104) 


D  (z)  -  1.5463(z  -  0.2934) 


If  Gc(s)  is  discretized  by  impulse  invariant  transformation,  then 


Dj/z) 


3z  -  1 
z  -  1 


(105) 


Case  2.  T  *  1  sec. 

Implementing  (76)  on  the  computer  program  gives 

.  -  1. 16415 (z  -  0.3679)  (z  +  0.3405)  (z  -  0.1362)  (z  +  4.3678)  nn,. 

(z  -  l)(z  -  0.2384)  (z  +  0.2978)  (z  +  4.0066)  UOt>j 

Discounting  the  pole  and  zero  outside  the  unit  circle  by  setting  z  =  1 
in  the  corresponding  factors  gives  the  stable  controller  as 


,  .  _  1.248(z  -  0. 3679) (z  +  0.3405Hz  -  0.1362) 
V }  (z  -  l)(z  -  0.2384) (z  +  0.2979) 


(107) 


The  prewarped  Tustin  transformation  and  impulse  invariant  transformation 
of  Gc(s)  at  T  *  1  are  given  by,  respectively. 


n  {z)  „  j.5574(z.,+  0.218) 
z  -  1 


(108) 


D^(z) 


3z  -  1 
z  -  1 


(109) 


Case  3.  T  *  2  sec. 


n(9s  ,  1.2145944(z  -  0. 1353) (z  +  0.2804) (z  -  0.0065) (z  +  4.7494)  ,1in. 
w  (z  -  l)(z  +  0.2701) (z  -  0.0221) (z  +  4.6059)  U  ' 


Discounting  the  unstable  pole  and  zero  gives: 


n,  *  .  1. 2656(z  -  0. 1353) (z  +  0.2804)(z  -  0.0065) 
K>  (z  -  l)(z  +  0.2701)(z  -0.0221) 


(111) 
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T 


The  prewarped  Tustin  transformation  and  impulse  invariant  transformation 
of  Gc(s)  at  T  ■  2  sec.  are,  respectively. 


Dt(z) 


Di(z) 


-1.185(z  +  2.687] 
(z  -  1) 

3z  -  1 
z  -  1 


(112) 

(113) 


The  comparison  of  unit-step  responses  of  the  closed-loop  systems 
using  D^(z),  Dt(z)  and  D(z)  are  shown  in  Fig.  6  and  Fig.  7.  The  com¬ 
parison  for  T  ■  2  is  not  plotted  because  the  performances  of  Dt(z)  and 
D^(z)  are  too  poor.  Fig.  8  compares  the  time  responses  of  the  state 
trajectory  matching  design  at  three  different  sampling  frequencies. 
Closed-loop  frequency  responses  are  shown  in  Figs.  9,  10  and  11. 

Example  2.  In  this  example  we  use  a  system  that  was  used  in  Reference 
II]  in  order  to  compare  results  with  the  frequency  response  matching 
method  for  digitalizing  the  continuous  controller.  Referring  again  to 
Fig.  1,  let  the  plant  and  the  continuous  controller  be 


10 

s(s  +  1) 


(114) 


1  +  24 
G  (s) - ^ 

C  i  +  -S_ 

1  7.2 


(115) 


respectively.  The  coefficients  in  the  state-variable  model  are: 


A  - 

L°  -±J 

b  ■ 

y 

(116) 

c  “ 

tl  0] 

d  - 

[0] 

(117) 

Ac  - 

[-7.2] 

*c  “ 

[7.2] 

(118) 

C  - 

c 

[-2] 

d  c 
c 

[3] 

(119) 
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Figure  6.  Comparison  of  Unit-Step  Responses  of  Example  1,  T 
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Figure  7 


Comparison  of  Unit-Step  Responses  of  Example  1,  T 


Output 


Figure  8.  Unit-Step  Responses  of  the  State-Trajectory-Matching  Design 
at  Different  Sampling  Frequencies 


We  choose 


and  consider  the  cases  of  T  ■  0.05,  T  *  0.15  and  T  =  0.30.  Note  that  for 
T  *  0.3,  the  sampling  frequency  is  only  slightly  above  twice  the  band¬ 
width  of  the  closed-loop  system.  For  this  system,  digitalization  by 
Tustin  transformation  gives  inferior  performances  both  in  the  time  domain 
and  the  frequency  domain,  and  can  be  found  in  [15] .  Hence,  comparison 
with  Tustin  method  will  not  be  presented.  Performance  comparisons  will 
be  made  with  the  frequency  response  matching  method  of  Reference  [1] . 

Case  1.  T  =  0.05 

Implementing  (76)  on  the  computer  program  gives 


2.8219(z  -  1) 


0. 9512) (z  +  0.2527) (z  -  0.8869) 


(z  -  0. 9995) (z  -  0. 9533) (z  +  0.2682)(z  -  0.9421) 


(z  -  0. 6742) (z  -  1.0514) (z  +  3.7310) 


(121) 


Simplifying  D(z)  by  cancelling  the  poles  and  zeros  that  are  very  close 
together  and  discounting  the  pole  and  zero  outside  the  unit  circle  (by 
setting  z  «  1  in  those  factors)  gives  the  reduced  version  of  D(z)  as 


2. 6898 (z  -  0. 9512) (z  +  0.2527)(z  -  0.8869) 
(z  -  0. 9421) (z  +  0. 2682) (z  -  0.6742) 


(122) 


Note  that  the  zero  at  0.9512  in  D(z)  is  also  a  pole  of  G.G(z),  the  z- 

1  -  e-Ts  h 

transform  of  -  G(s). 

8 

An  attempt  has  been  made  to  obtain  a  digital  controller  via  frequency 
response  matching  method  of  Reference  [1]  for  comparison.  But  the 
method  failed  to  provide  a  stable  digital  control  system  for  T  *  0.05. 

Case_2 .  T  ■  0. 15 

Implementing  (76)  on  the  computer  program  for  T  -  0.15  yields 
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D(z) 


2.6577(z  +  0.2265)  (z  -  0.6977)(z  -  0.8607)2 
(z  -  0. 2040) (z  +  0.2771) (z  -  0. 8607) (z  -  0. 8090) 


(z  -  l)(z  -  1. 1618) (z  +  3.0744) 
X  (z  -  l)(z  -  1.1618) (z  +  3.7062) 


(123) 


The  reduced  model  of  D{z)  is 


D(z) 


2.294(z  +  0. 2265)  (z  -  0.6977Hz  -  0.8607) 
(z  +  0.2771) (z  -  0. 2040) (z  -  0.8090) 


(124) 


In  applying  the  frequency  response  matching  method,  one  can  specify  the 
order  of  the  digital  controller  desired.  The  first  and  second  order 
digital  controller  obtained  by  the  frequency  response  matching  method 
are  [1], 


(  \  3. 6787(z  -  0.6275) 

DflU'  "  z  +  0.3837 


and 


Df2(Z) 


3.96  -  3. 65z-1  +  0. 56z 
1  +  0. 34z-1  -  0. 384z~ 


(125) 


(126) 


Results  of  similation  shows  that  while  the  frequency  domain  performance  of 
the  second  order  controller  seems  to  be  slightly  better  than  the  first 
order  controller,  ratcheting  and  jittering  begin  to  appear  in  the  time 
response  of  the  system  using  the  second  order  controller.  Hence  in 
order  not  to  overcrowd  the  graphs,  the  second  order  controller  is  regarded 
as  undesirable  and  not  used  in  the  performance  comparison,  and  higher 
order  controllers  by  the  frequency  response  matching  method  are  not 
sought . 

Case  3.  T  ■  0.30 

The  transfer  function  of  the  digital  controller  (Eq.  (76))  is  found 

to  be 
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2.0175(z  +  0. 1958)  (z  -  0.4865)(z  -  0.7394) 
(z  +  0.1099) (z  +  0.3678) (z  -  0.6012) 


D(z)  = 


(z  -  0. 7423) (z  -  l)(z  -  1. 3499) (z  +  2.4504) 
X  (z  -  0. 7414) (z  -  l)(z  -  1.3499) (z  +  3.5534) 


(127) 


The  reduced  model  of  D(z)  is 


D(z) 


1.5288(z  +  0.1958)  (z  -  0.4865)(z  -  0.7394) 
(z  +  0. 1099) (z  +  0. 3678) (z  -  0.6012) 


(128) 


For  T  =  0.30,  the  frequency  response  method  fails  to  produce  stable 
digital  controller  again. 

The  unit-step  responses  at  the  output  of  the  state  trajectory 
matching  design,  the  continuous  model  and  the  system  obtained  by  fre¬ 
quency  response  matching  method  (first-order  controller  only)  for  T  = 

0.15  are  compared  in  Fig.  12.  Fig.  13  compares  the  squared  errors  of 
the  time-responses  of  the  state  trajectory  matching  design  with  the  fre¬ 
quency  response  matching  design  using  the  continuous  model  as  the  ideal 
response.  The  time  responses  of  the  second  state  variable  y  of  the  two 
designs  at  T  *  0.15  are  shown  in  Fig.  14.  The  squared  errors  are  shown 
in  Fig.  15.  Comparisons  of  time  domain  performances  of  the  state  trajec¬ 
tory  matching  design  for  different  sampling  frequencies  are  made  in 
Fig.  16  and  Fig.  17.  The  frequency  response  comparisons  are  made  in 
Fig.  18  and  Fig.  19. 

It  should  be  noted  that,  although  a  first  glance  of  Fig.  12  may 
give  the  impression  that  the  frequency  matching  method  and  the  state- 
trajectory  matching  method  both  give  about  the  same  time-domain  per¬ 
formances,  Fig.  13  shows  that  the  state-trajectory  matching  gives  a 
much  better  design.  Fig.  14  and  Fig.  15  show  that  the  state  matching 
method  yields  much  better  step  response  in  the  other  state  variable  too. 
Although  the  frequency  matching  method  seems  to  yield  a  frequency  re¬ 
sponse  that  matches  the  frequency  response  of  the  continuous  model  better 
at  the  higher  frequency  range  than  that  of  the  state  matching  design 
(Fig.  19),  Appendix  D  suggests  that  the  close  matching  of  frequency  re¬ 
sponses  between  a  discrete  time  system  and  a  continuous  system  at  the 


36 


higher  frequency  range  is  not  only  meaningless,  but  also  undesirable  if 
the  discrete  time  system  is  actually  a  digitally  controlled  continuous 
system. 


i 
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Figure  12.  Comparison  of  Unit-Step  Responses  of  Example  2  for  T 
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Figure  13.  Squared  Errors  (Relative  to  the  Continuous  Model)  of  the 
Time  Responses  of  Figure  12 
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Figure  14.  Unit-Step  Responses  of  the  Second  State  Variable  of 
Example  2  at  T  ■  0.15 
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Figure  15.  Squared  Errors  (Relative  to  the  Continuous  Model)  of  the 
Time  Responses  of  Figure  14 
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Figure  17.  Unic-Step  Responses  of  Che  Second  SCate  Variable  y  of  Che 

Scace-Trajeccory-Macching  Design  of  Example  2  for  DifferenC 
Sampling  Frequencies 
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Frequency  Response  Comparison  for  Two  Different  Designs  of  Example  2  fi 
Sampling  Period  T  =  0. 15  sec. 


SECTION  VI 


VERIFICATION  OF  THE  OPTIMAL  CONTROL  SEQUENCE 
VIA  THE  DISCRETE  MAXIMUM  PRINCIPLE 

It  has  been  shown  in  Section  III  that  the  extended  maximum  prin¬ 
ciple  offers  an  efficient  way  to  derive  the  optimal  control  sequence  for 
a  sampled-data  system.  The  extended  maximum  principle,  though  powerful 
in  this  and  many  other  applications  [5-9]  has  not  been  widely  used  in  the 
technical  community.  In  many  cases  where  the  extended  maximum  principle 
may  be  used,  an  equivalent  discrete  performance  index  is  derived  and 
the  discrete  maximum  principle  is  used  instead  [12-14], 

It  will  be  shown  in  this  Section  that  the  optimal  control  sequence 
of  Eq.  (42)  may  be  derived  by  using  the  discrete  maximum  principle  but 
a  more  cumbersome  procedure  is  involved. 

The  performance  index  (5)  may  be  written  as 
.  N-l  (k+l)T 

J  =  £2  i  1  5  t*'  (t)Qx(t)  -  2x’(t)Qx  (t)  +  x*(t)Qx  (t) 

N  2  k=0  kT  ^  a 

+  Bu2(t)  -  2Bu(t)u  (t)  +  0u2(t)]dt  (129) 

m  m 

Since  u(t)  is  the  output  of  a  zero-order-hold,  substituting  (10),  (28), 
(30)  into  (129)  (except  for  the  term  x^(t)Qxa(t))  gives: 

14_  i  N-l 

J  “  N^»  2  1  lQ(T)2E.0tT>  +  2M(T)bu(kT) 

-  2Qm(T)xm(kT)  -  2Mm(T)bma]  +  6(T)u2(kT) 

(k+l)T 

-  2u(kT) [b'M  (T)x  (kT)  +  6m(T)a  +6  /  um (t)dt] 

~  m  — m  m  ^  m 

(k+l)T  (k+l)T  . 

+  J  x  (t)Qx  (t)  +  0  /  u  (t)dt)  (130) 

kT  —a  —a  kT  m 
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where 


Q(T) 

M(T) 

6(T) 

VT) 

Mm(T) 

Mm(T) 

Bm(T) 


T 

/ 

0 

T 

/ 

0 

8T 

T 

/ 

0 

T 

/ 

0 

T 

/ 

0 


b’ 


$' (t)Q$(t)dt 

♦ ' (t)Q$S(t)dt 
T  , 

+  b '  /  $S  (t)Q$S(t)dt 
0 

*’(t)Q  $®(t)dt 
m  m 

«S’(t)Qm«m(t)dt 

T 

/  *S'(t)Q  *®(t)dtbn 

q  mm  -in 


(131) 

(132) 

(133) 

(134) 

(135) 

(136) 

(137) 


The  discrete  time  system  to  be  optimized  is  described  by  (47),  which  is 
repeated  here  for  convenience: 


x[(k  +  1)T]  -  $(T)x(kT)  +  $s(T)bu(kT) 


(138) 


In  the  formulation  of  discrete  maximum  principle,  the  Hamiltonian  is 
given  by  {12]  (for  the  sake  of  notational  brevity,  the  sampling  period 
T  is  dropped  in  the  argument  of  the  discrete-time  functions): 


H(k)  -  H(x(k) ,  £(k  +  1),  u(k)) 

-  £'  (k  +  l)[$(T)x(k)  +  $s(T)bu(k) ] 

+  j  x’(k)[Q(T)x(k)  +  2H(T)bu(k)  -  2Qm(T)xm(k) 

-  2M  (T)b  a]  +  ^  6(T)u2(k)  -  u(k) [b'M(T)xm(k) 

m  —m  l  tn  tu 
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(k+l)T  .  (k+l)T 

+  Bn(T)a  +  8  /  +jf  xa(t)Qxa(t)dt 

kX  kT 

1  /m)T  2 

+  f  8  /  u^(t)dt  (139) 

kT 


The  difference  equation  for  the  adjoint  vector  is 


£<k) 


3H(k) 

3x(k) 


-  Q(T)x(k)  +  M(T)bu(k)  -  (^x^k) 

-  M  (T)b  a  +  #' (T)p(k  +  1) 

m  — m 


(UO) 


The  coupling  equation,  or  the  condition  of  optimality  is  given  by 


3H(k) 

3u(k) 


0 


(141) 


Substituting  (139)  into  (141)  gives 

£'  (k  +  l)$S(T)b  +  x* (k)M(T)b  +  B(T)u(k) 

(k+l)T 

-  b'M  (T)x  (T)  -  B  (T)o  -  B  /  u  (t)dt  -  0  (142) 

kT 

Solving  the  above  equation  for  u(k)  yields  the  optimal  control  sequence: 


u(k)  -  (b'M  (T)x  (k)  -  b'M' (T)x(k)  +  8  (T) a 

„  /  —  m  -m  —  —  m 


B(T) 


«•  (kfl)T 

-  b'«  (T)£(k  +  1)  +  B  J  u  (t)dt] 

kT  m 


(143) 


If  it  is  recognized  that 


*'(’t)  -  eA’T  -  ip(-T)  -  [^(T)]”1 


(144) 
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then  p(k  +  1)  may  be  obtained  from  (140)  as 


p(k  +  1)  -  <|»(T)[£(k)  -  Q(T)x(k)  -  M(T)bu(k) 

+  C^COx^k)  +  Mm(T)bma]  (145) 

Also,  in  view  of  Eqs.  (22),  (29)  and  (31), 

(k+l)T  _ 

/  u  (t)dt  *  c  $®(T)x  (k)  +  c  $SS(T)b  a  +  d  aT  (146) 

m  mo  — m  m  m  - m  m 

Substituting  (145)  and  (146)  into  (143)  and  solving  for  u(k)  give  the 
explicit  expression  of  the  optimal  control  sequence  u(k)  as 

u(k)  -  - - -r -  {-b’4S,(T)iMT)£(k) 

B(T)  -  b'$  (T)i^(T)M(T)b 

+  [-b'M’(T)  +  b'*8’ (T)4»(T)Q(T)]x(k) 

+  rb'Mm(T)  +  0cm*8(r)  -  b,*S'(T)^(T)Qm(r)]2Em(k) 

+  (6  (T)  +  0c  «SS(T)b  +  Bd  T  -  b’*8'  (T)i|»(T)M  (T)b  ]a}  (147) 

TO  TO  121  “TO  TO  —  TO  "“TO 

It  remains  to  be  shown  that  the  right-hand  side  of  (147)  is  equal  to  the 
right-hand  side  of  (42).  First,  some  preparation  is  in  order.  It  is 
seen  that,  in  view  of  (132)  and  (144) 

T  t 

M'  (T)  -  /  /  iK-A)dAQ*(T)dx 
0  0 

T  T 

■  /  /  ♦(*  “  x)dxQ*(x)dx 

0  0 

T  T 

*  /  /  'Kx  “  x)Q$(x)dxdx  (148) 

0  x 


49 


Note  that  a  change  of  variables  and  a  change  of  order  of  integrations  are 
involved  in  obtaining  the  last  expression  in  (148).  Furthermore,  it  is 
seen  from  (38)  and  (39)  that 


FS(T)  ■  /  /  i|/(x  -  T)Q$(x)dxdx 

0  0 


(149) 


Thus 

T  T 

F  (T)  +  M'(T)  -  /  /  i l>(x  -  x)Q$(T)dTdx 

0  0 

T  T 

■  /  ij»(x  -  T)J  i|i(T  -  x)Q*(x)dxdx 
0  0 


T 

-  /  i|«(-X)F(T)dX 
0 

-  #S'(T)F(T)  (150) 


Therefore 

FS(T)  -  -M’ (T)  +  48' (T)F(T)  (151) 


Now  changing  the  order  of  integration  in  the  first  expression  of  M'  (T)  in 
(148),  with  the  aid  of  (144)  and  (132),  gives 


T  x 

M'(T)  -  /  /  iK-X)dXQ$(x)dx 
0  0 


T  T 

"  /  /  ’j'(“X)Q«(x)dTdX 

0  X 


T  T 

■  /  ♦(-X)Q /  #(x)dxdX 
0  X 


T  T 

-  /  i|* (— X)Q  /  ♦(x)dxdX 
0  0 


T  X 

-  /  <j»(-X)Q  /  ®(x)dTdX 
0  0 
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(152) 


-  *S'(T)Q*S(T)  -  /  i|»(-X)Q*8(X)dX 

0 

-  4>S'(T)Q$S(T)  -  M(T) 


Substituting  (152)  Into  (151)  gives 

FS(T)  -  -4S'(T)Q*S(T)  +  $S' (T)F(T)  +  M(T) 


(153) 


But  since  <j;(T)  is  a  matrix  exponential,  M(T)  may  be  expressed  In  the 
following  form,  with  the  aid  of  the  theorem  given  in  Appendix  C: 


T 

M(T)  -  /  <t'(t)Q*S(t)dt 
0 

T  t 

-  *(-T)  /  <fr(T  -  t)Q  /  «(X)dXdt 
0  0 

T  t 

“  «K-T)  /  /  ^ (t  -  T)Q*(r)dTdt 

0  0 

T 

.  ^(-t)  /  F(t)dt  -  iK-T)FS(T)  (154) 

0 


Hence,  in  view  of  (144)  and  (154),  FS(T)  of  (153)  may  be  written  as: 

FS(T)  ■  -$S(T)Q*S(T)  +  $8’ (T)F(T)  +  (T)FS(T)  (155) 


Integrating  both  sides  of  (155)  from  0  to  I  gives 


F88(T)  «  -  /  *8'(t)Q«8(t)dt  +  *8’(T)F8(T) 
0 


(156) 


Now  it  can  be  shown  that  the  control  sequence  given  in  (147)  is 
identical  to  that  given  in  (42).  With  the  aid  of  (133),  (154)  and  (156)  the 
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denominator  in  the  right-hand  side  of  (147)  may  be  reduced  as  follows: 


B(T)  -  b’*S’(T)<|t(T)M(T)b 
T 

-  ST  +  b’ [/  *®'(t)Q*®(t)dt  -  *S'(T)iKT)M(T)]b 

0 

T  ,  , 

-  ST  +  b' [/  «8  (t)Q$S(t)dt  -  *S  (T)FS(T)]b 

0 

-  ST  -  b’FSS(T)b  (157) 

Because  i|/(T)  is  a  matrix  exponential,  the  coefficient  of  j>(kT)  in  (147)  may 
be  reduced  as  follows: 

,  T  T 

«  (T)iJi(T)  -  /  d»(-t)dtip(T)  -  /  <K T  -  t)dt 
0  0 

T 

-  /  i|)U)dX  -  At)  (158) 

0 

The  coefficient  of  x(k)  may  be  reduced  as  follows: 

-M'(T)  +  4S'(T)<KT)Q(T) 

T 

-  -M’(T)  +  48'  (T)ifr(T)  /  *(- t)Q4(t)dt 

0 

-  -M’(T)  +  4®'(T)F(T)  -  FS(T)  (159) 

where  (131),  (38)  and  (151)  have  been  used.  Similarly  the  coefficient  of 
x  (k)  in  (147)  is  shown  to  be 

“ID 

b'Mn(T)  +  6cm«8(T)  -  b,48'(T)<i(T)Qm(T) 

-  Sc  At)  -  b'F®(T)  (160) 

tp  m  id 
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For  Che  coefficient  of  a  in  (147),  it  is  seen  that 


6  (T)  -  b'«S'(T)i|»(T)M  (T)b 

m  m  "tii 

T 

-  b '[/  •8'(t)<^*“(t)dt  -  *S'(T)F8(T))bm 

-  -b'FSS(T)b  (161) 

—  id  -m 

ss 

Where  formulas  (135),  (40),  (137)  and  an  equivalent  of  (156)  for  F^  (T) 
have  been  used. 

Finally,  the  following  relations  (which  have  been  found  in  the  deriva¬ 
tions  of  (154)  -  (161)) 


<J>(T)Q(T)  -  F(T)  (162) 

*(T)M(T)  -  FS(T)  (163) 

<KT)Qo(T)  -  Fffl(T)  (164) 

i|)(T)M  (T)  -  FS(T)  (165) 

m  m 


lead  to  the  fact  that  the  solution  £(k)  of  (145)  is  the  negative  of  the 
solution  of  the  difference  equation  given  by  (49).  This  fact  and  Eqs. 

(157)  -  (161)  imply  that  the  optimal  control  sequence  given  in  (147)  by 
the  discrete  maximum  principle  is  identical  to  the  optimal  control  sequence 
given  in  (42)  by  the  extended  maximum  principle.  It  is  also  clear  that 
the  derivation  of  (42)  is  much  simpler  than  the  derivation  of  (147).  The 
extended  maximum  principle  even  yields  more  concise  expressions. 
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SECTION  VII 


CONCLUSIONS  AND  RECOMMENDATIONS 


The  extended  maximum  principle  has  been  applied  to  the  problem  of 
designing  a  digital  controller  that  drives  the  state  trajectory  of  a 
continuous  plant  in  an  output  feedback  system  to  follow  its  continuous 
model  as  closely  as  possible.  Closed  form  expressions  of  the  optimal 
digital  controller  have  been  obtained,  as  given  by  Equations  (74),  (76), 
(90)  and  (95).  The  result  has  been  verified  by  a  derivation  via  the 
widely  used  discrete  maximum  principle,  which  would  entail  a  far  more 
cumbersome  procedure  to  obtain  the  same  result. 

The  optimal  digital  controller  gives  rise  to  a  closed-loop  system 
(Fig.  2  or  Fig.  4)  that  has  a  discrete  transfer  function  equal  to 
(Eq.  (82)  repeated  here): 


21*1 

R(z) 


b'H(z)b 

m 

b’H(z)b 


GhG(z) 


(166) 


if  8  -  0  (integral  squared  error  of  the  control  signal  is  not  weighted  in 
the  performance  index).  A  similar  expression  may  be  derived  for  6  J*  0. 
The  formulas  for  the  digital  controller,  given  by  (74),  (76),  (90)  and 
(95),  are  rather  difficult  to  compute.  Numerical  computations  of  these 
equations  always  result  in  poles  and  zeros  in  the  unstable  region  that 
should  not  be  implemented.  Owing  to  the  amount  of  computations  involved, 
it  Is  not  clear  whether  the  unstable  poles  and  zeros  in  the  digital  con¬ 
troller  should  be  identical  to  one  another  and  thus  be  cancelled,  or  are 
introduced  by  z- transforming  the  control  and  error  sequences.  More  study 
is  needed  in  this  direction. 

For  the  lower  order  systems  (second  and  third  order  systems),  this 
method  of  digital  controller  design  yields  excellent  results,  even  though 
the  sampling  frequencies  are  relatively  low.  The  digital  controllers 
designed  by  the  present  method  compares  favorably  with  those  designed  by 
the  frequency  matching  method,  and  are  far  superior  to  those  obtained  by 
the  conventional  Tustln  or  impulse  invariant  transformations. 
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Attempts  have  also  been  made  to  compute  the  digital  controller  of 
Eq.  (95)  for  a  sixth  order  model  of  the  Inner  loop  of  the  longitudinal 
control  of  YF-16  fighter  aircraft  at  30,000  ft  altitude  and  Mach  0.6. 

It  is  found  that  (95)  is  not  well  conditioned  for  numerical  computation 
for  higher  order  systems  in  the  sense  that  it  involves  ratios  of  differ¬ 
ences  between  polynomials  which  are  approximately  equal  to  one  another. 

This  and  the  fact  that  a  large  amount  of  multiplication  and  integration 

are  required  in  computing  H (z)  and  H  (z)  make  it  impossible  to  obtain  a 

m 

valid  expression  of  the  digital  controller  by  way  of  the  computer  programs 
available  to  or  developed  in  this  project.  Another  difficulty  involved 
in  computing  (95)  is  the  high  order  in  the  resulting  transfer  function. 

The  ill-conditionedness  of  the  digital  controller  formulas  may  perhaps 
be  circumvented  by  transforming  the  z-plane  formulas  into  the  w' -plane,  i.e. , 


z 


1  +  w 


1  -  w 


(167) 


The  w'-plane  concept  was  initiated  by  Whitbeck  [16].  W' -transformation 
has  the  advantage  over  the  w  transformation  in  that  when  T  -*■  0,  w'  -*•  s. 

It  may  also  be  easier  to  find  reduced-order  models  of  the  digital  con¬ 
troller  in  the  w'-plane  in  order  for  this  method  to  apply  to  practical 
higher  order  systems. 

Equation  (166)  represent  the  ideal  closed-loop  transfer  function 
when  the  optimal  digital  controller  is  employed.  It  represents  the  best 
closed-loop  digital  control  system  that  can  be  achieved  if  the  digital 
system  is  required  to  model  after  the  continuous  system. 

The  numerical  computation  of  the  optimal  digital  controller  transfer 
function  of  Eqs.  (74),  (76),  (90)  and  (95)  involves  integrations  and 
double  integrations  of  matrix  exponentials  and  of  convolutions  of  matrix 
exponentials.  We  used  Romberg's  algorithm  [17]  on  an  IBM  370  computer 
with  extended  precision.  For  a  sixth  order  system,  the  CPU  time  required 
to  compute  a  set  of  coefficients  for  (76)  is  0.25646  hours  for  T  ■  0.04, 
and  0.74129  hours  for  T  -  0. 1.  Therefore,  a  more  efficient  computer  program 
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X  d£ 

such  as  one  that  performs  analytic  integration  by  giving  —  e  as  the 

at  ^ 

integration  of  e  for  matrix  functions  needs  to  be  developed.  The 

analytic  integration  program  would  not  only  cut  down  the  CPU  time,  but 
also  improves  precision  since  numerical  values  are  substituted  at  the  last 
stage.  It  is  useful  not  only  for  computing  the  formulas  of  the  digital 
controller  developed  in  this  project,  but  also  in  discretizing  a  contin¬ 
uous  control  system  or  in  applying  the  discrete  maximum  principle  to 
sampled-data  systems. 
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APPENDIX  A 


THE  OVERALL  STATE  EQUATIONS  OF  THE  CONTINUOUS  MODEL 


The  feedback  structure  in  Fig.  1  gives 


e  (t)  *■  r (t)  -  y  (t) 
id  m 


(A-l) 


Eliminating  e  (t)  and  y  (t)  among  Eqs.  (A-l),  (2)  and  (4)  and  solving 
m  m 

for  u  (t)  gives 
m 


-d  c 

u  (t) - - 

m 


l+dd  ^a(t)  +  1+dd  +  1+dd  r(t) 

c  c  c 


(A-2) 


Eliminating  u  (t)  and  y  (t)  among  Eqs.  (A-l),  (2)  and  (4)  and  solving 
m  in 

for  e  (t)  gives 
m 


em(t)  "  1+dd 


-dc  . 

+  1+dd  +  1+dd  r(t) 

c  c  c 


(A-3) 


Substituting  (A-2)  into  (1)  and  (A-3)  into  (3)  yields,  respectively. 


.<*> 


/r  bd  c  \ 


be 

2^(0  +  C 


bd 


1+dd  ^c(t)  +  1+dd  r(t) 
c  ** 


(A-4) 


i  (t) 


-b  c 
—c 

1+dd 


x  (t)  +  (  A  - 


(b  dc  \ 

AC  ■  3dd^) 


X  (t)  + 


b 

— c 


1+dd 


r(t) 


(A-5) 
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APPENDIX  B 


SOLUTIONS  OF  THE  STATE  AND  ADJOINT  EQUATIONS 


The  solution  £(t)  of  Eq.  (14),  in  terms  of  x(t)  and  x  (t),  is 


given  by 


£(t)  =  i|«(t-kT)£(kT)  +  /  \j)(t-T)Q[x(t)  -  x  (x)]dx 

kT 


(B-l) 


for  kT  <  t  <  kT,  where 


\Ji(t-kT)  =  e 


-A’ (t-kT) 


(B-2) 


With  the  aid  of  Eq.  (28),  Eq.  (B-l)  can  be  written  as 


£(t)  ■  ij»(t-kT)£(kT)  +  /  ip  (t— t )Qx(t )dx 

kT 


-  /  tKt-T)QnXm(T)dr 


(B-3) 


The  vectors  x  (t)  and  x(t)  are  solutions  of  Eq.  (21)  and  Eq.  (7), 
— m  — 

respectively.  They  are  given  by,  for  kT  <  t  <  (k+l)T, 


x  (t)  -  *  (t-kT)x  (kT)  +  $  (t-kT)b  a 
— m  m  “in  m  m 


x(t)  -  $ (t-kT)x(kT)  +  <r (t-kT)bu(kT) 


(B-4) 

(B-5) 


where 


A  (t-kT) 


»  (t-kT)  -  e 

IQ 


m 


(B-6) 


t  t-kT 

<&S(t-kT)  ■  /  $  (t-x)dx  ■  /  $  (x)dx 

jjj  •  m  *  m 


ki  m  o 


(B-7) 


$(t-kT)  -  eA(t_kT) 


(B-8) 
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t  t-kT 

*  (t-kT)  -  /  $(t-x)dx  -  /  ® (x)dx 

kT  0 


(B-9) 


Now  substituting  (B-4)  and  (B-5)  into  (B-3)  and  invoking  (B— 7)  and 
(B-9)  gives,  for  kT  <  t  <_  (k+l)T, 


£(t)  =  i|.(t-kT)£(kT)  +  /  <Kt-T)Q4(t-kT)dTx(kT) 

kT 

t  X 

+  /  <|i(t-x)Q  /  <Kx-X)dXdxbu(kT) 

kT  kT 


-  /  ij/ ( t— r) Q_4>  (x-kT)dxx  (kT) 

,  L  in  in  —m 


kT 


t  T 

-  /  iKt-r)Q  /  $  (x-X)dXdxb  a 

'  m  ,  '  in  — m 

kT  kT 


(B-10) 


In  view  of  the  identities  derived  in  Appendix  C,  a  new  set  of  notations 
may  be  used  for  the  integrals  which  appear  in  the  above  equation,  for 
the  sake  of  clarity: 

£(t)  =  i|>(t-kT)£(kT)  +  F(t-kT)x(kT)  +  FS(t-kT)bu(kT) 


F  (t-kT)x  (kT)  -  Fs(t-kT)b  a 
m  — m  m  — tq 


(B-ll) 


where 


F(t-kT)  -  /  ^(t-x)Q$(x-kT)dx 
kT 

t-kT 

-  /  i|»(t-kT-T)Q«(x)dx 
0 


(B-12) 
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t  T 

FS(t-kT)  -  /  iKt-x)Q  /  $(x-X)dXdx 

kT  kT 

t-kT  x 

-  /  /  *(x-X)Q$(X)dXdx 

0  0 

t-kT 

-  /  F(t)dT 


F  (t-kT)  =  /  i|»(t-T)Qm$!n(T-kT)dT 

m  kT 

t-kT 

-  /  iKt-kT-T)<^*m(x)dx 


F  (t-kT) 
m 


t  T 

/  <Kt— t)Q  /  *m(t-X)dXdT 

kT  kT 

t-kT  t 

/  /  ♦(x-X)Q  •  (X)dXdx 

0  0 


t-kT 

/  F  (x)dx 
0 


(B-13) 


(B-14) 


(B-15) 


FORMULAS  OF  FINITE  MULTIPLE  INTEGRALS 


THEOREM;  The  convolution  between  f(t)  and  the  area  under  g(t),  if  it 
exists,  is  equal  to  the  area  under  the  convolution  between  f(t)  and 
g(t). 

T  x  T  t 

/  f(T-x)  /  g(X)dXdx  8  |  /  f (t-x)g(x)dxdt  (C-l) 

0  0  0  0 

Proof. 


In  view  of  Fig.  C-l,  interchanging  the  order  of  integration  for  the 
integral  in  the  right-hand  side  of  Eq.  (C-l)  gives 


T  t  T  T 

/  J  f (t-x)g(x)dxdt  -  /  /  f (t-x)g(x)dtdx 


(C-2) 


Substituting  t*  for  t-T  and  dt*  for  dt  (regarding  r  as  a  constant  in  the 
partial  integration  with  respect  to  t)  in  Eq.  (C-2),  we  obtain,  after 
changing  the  integration  limits  accordingly. 


T  T  T  T-t 

/  /  f (t-x)g(x)dtdx  -  /  /  f (t')dt'g(T)dT  (C-3) 

0  t  0  0 


In  a  similar  fashion,  interchanging  the  order  of  integration  in  the 
left-hand  side  of  Eq.  (C-l)  gives 


T  x  T  T 

/  f (T-t)  /  g(X)dXdx  -  /  /  f (T-T)g(X)dxdX  (C-4) 

0  0  OX 


Substituting  t  for  T-t  and  dt  for  -dr  and  changing  the  integration  limits 
accordingly  in  the  right-hand  side  of  Eq.  (C-4)  gives 


T  T  TO 

/  /  f(T-T)g(X)dxdX  -  /  J  (-l)f(t)dtg(X)dX  (C-5) 

OX  0  T-X 


Comparison  of  the  Eq.  (C-3)  with  Eq.  (C-5)  completes  the  proof. 

The  following  is  a  list  of  formulas  that  are  used  in  the  deriva¬ 
tions  in  Appendix  B.  The  proof  involves  changes  of  variables  or 
changes  of  order  of  integrations,  and  is  hence  omitted. 


Formula  1. 


Formula  2. 


Formula  3. 


Formula  4. 


Formula  5. 


t  t-kT 

/  f(t-x)dx  -  /  f(x)dx 
kT  0 

t  t-kT 

/  f (t-x)g(x-kT)dx  -  /  f (t-kT-x)g(T)dx 

kT  0 

t  x  t-kT  x 

/  f(t-x)  /  g(x-X)dXdx  -  /  f(t-kT-T)  /  g(X)dXdx 
kT  kT  0  0 

t  x  t-kT  x 

/  f(t-x)  /  g(T-X)dXdx  -  /  /  f (x-X)g(X)dXdx 

kT  kT  0  0 

(k+l)T  t  T  t 

/  /  f (t-x)g(x-kT)dxdt  -  /  /  f (t-x)g(x)dxdt 

kT  kT  0  0 


64 


Formula  6. 


Formula 


Formula 


Formula 


Formula 


t  t  t  x 

/  f(t-x)/  g(x-A)dAdx  -  /  f(t-T )/  g(A)dAdx 

0  0  0  0 

t  X  t  T 

7.  /  f(t-x)/  g(A)dAdx  -  /  /  f(x-A)g(A)dAdx 

0  0  0  0 

t  x  t  x 

8.  /  f(t-x)/  g(x-X)dXdx  ■  /  /  f (x-X)g(X)dXdx 

0  0  0  0 

(k+l)T  t  x 

9.  /  /  f(t-x)  /  g(x-X)dXdxdt 

kT  kT  kT 

T  t  x 

•ff  f(t-x)/  g(X)dXdxdt 
0  0  0 

(k+l)T  t  x 

10.  /  /  f(t-x)  /  g(x-X)dXdxdt 

kT  kT  kT 

T  t  x 

”  /  /  /  f<T-X)g(X)dXdxdt 

0  0  0 
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APPENDIX  D 

SUBALIASES  IN  THE  FREQUENCY  RESPONSE 
OF  DIGITALLY  CONTROLLED  SOURCES 


SUBALIASES  IN  THE  FREQUENCY  RESPONSE 
OF  DIGITALLY  CONTROLLED  SYSTEMS 


Hsl-Han  Yeh* 

University  of  Kentucky,  Lexington,  Ky. 
and 

Richard  F.  Whitbeck** 

Systems  Technology,  Inc.,  Hawthorne,  Calif. 


I.  INTRODUCTION 

In  a  recent  paper,  Whitbeck,  Didaleusky  and  Hofmann  [1]  extended 
the  concept  of  the  tradiational  "sampled  spectrum"  frequency  response 
for  discretely  excited  continuous  systems.  When  a  sinusoidal  wave  is 
input  to  a  discretely  excited  continuous  system,  N  sine  waves  at  dif¬ 
ferent  alias  frequencies  are  required  to  match  the  continuous  steady- 
state  time  response  at  the  sampling  instants  and  at  N-l  equally  spaced 
intersample  points.  In  the  special  case  of  N*l,  this  reduces  to  the 
traditional  concept  of  the  sampled-spectrum  frequency  response  for 
sampled-data  systems.  Letting  N  approach  infinity  gives  an  infinite 
spectrum  for  the  continuous  steady-state  response  of  a  discretely 
excited  continuous  system.  This  theory  enables  one  to  write  an  exact 
expression  of  the  time  response  sampled  at  any  rate  that  is  an  integer 
multiple  of  the  sampling  rate  of  the  system.  The  practical  value  of 
knowing  such  an  expression  of  the  output  is  evident. 


^Associate  Professor,  Department  of  Electrical  Engineering 
^Principal  Research  Engineer 


However,  in  the  derivation  of  Whitbeck,  Didaleusky  and  Hofmann 
[1],  only  positive  aliases  of  the  input  frequency  are  included  in  the 
representation  of  the  sampled  continuous  output.  Whereas  the  result  is 
correct  for  the  case  where  N  is  finite,  the  limit  case  where  N 
approaches  infinity  is  in  error  and  it  contains  only  half  of  the 
spectral  components  necessary  for  an  asymptotic  representation  of  the 
continuous  output.  In  this  paper,  a  new  derivation  which  includes 
subaliases  in  the  spectral  representation  of  the  sampled  output  will  be 
presented.  It  will  be  shown  that,  as  the  output  sampling  rate  N 
approaches  infinity,  the  infinite  spectrum  of  the  continuous  output 
contains  all  aliases  and  subaliases  of  the  input  frequency.  Also  pre¬ 
sented  will  be  a  direct  derivation  of  the  infinite  spectrum  of  the 
continuous  output,  without  invoking  the  expression  of  the  sampled 

steady-state  response.  This  confirms  the  correct  representation  of  the 

% 

continuous  steady-state  output  of  discretely  excited  continuous  systems 
in  response  to  a  single  sinusoidal  input.  This  derivation  brings  forth 
a  unified  concept  of  frequency  response  which  enables  one  to  write  the 
spectral  representation  of  the  output  of  a  discretely  excited  continuous 
system  on  the  basis  of  the  frequency  response  of  the  continuous  system. 

II.  FREQUENCY  COMPONENTS  IN  THE  SAMPLED  OUTPUT 

Consider  the  system  of  Fig.  1  where  G(s)  represents  an  arbitrary 
transfer  function  and  M(s)  represents  an  arbitrary  data  hold.  Let  the 
input  be  a  unit  amplitude  exponential  and  the  output  be  sampled 

with  period  T/N.  Using  multirate  sampling  results  (See  Appendix  of 
[1])  yields 


(1) 


.T/N 


[gmrt] 


T/N 


(GM)T/N  rt 


-  (GM) 


T/N  _ z 


N 


N  jbT 
z  -  e 


,  A  sT/Nv 
(z-  e  ) 


where  the  nocation  follows  Whitbeck,  Didaleusky  and  Hofmann  [1].  The 

superscript  denotes  the  period  of  sampling  operation.  In  the  time 

T  T 

domain,  a  sampled  function  [r(t)]  or  r  is  defined  by 


[r(t)]T  £  l  r (n) S(t-nT)  -  r(t)  \  6(t-nT) 
n«-«  n->-“ 


(2) 


In  the  s-domain,  a  sampled  function  R  is  defined  by 


„T,  .  A  c  ,  .  -nsx  x  r  ,  .xirn. 
R  (s)  “  I  r(n)e  m  j  l  R(s  +  j-y) 

n=0  n*-*  1 


-nsT  1 


,27m. 


(3) 


and  in  z-domain. 


_T.  .  A  _T,  . 
R  (z)  =  R  (s) 


sT 

e  -+-z 


(4) 


But  in  the  presence  of  a  higher  rate  sampler  with  period  T/N,  as  in 

_rp  y 

this  paper,  it  is  desirable  to  use  z-e  and 


RT(zN)  £  RT(s) 


sT/N 
e  -*z 


(5) 


Hence  it  is  simpler  to  think  of  frequency  response  of  discretely  excited 
continuous  systems  in  terms  of  s  rather  than  z.  Where  no  confusion  may 
arise,  the  variables  a  and  z  may  be  omitted  to  give  versatility  to  the 
notation.  For  the  sake  of  notational  brevity,  in  the  following  devel¬ 
opment  G(s)M(s)  will  be  written  as  GM(s),  for  occasions  where  a  needs  to 
be  substituted  by  a  string  of  notations. 
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The  steady-state  component  in  the  sampled  output  is  of  interest. 
In  taking  the  partial  fraction  expansion  of  the  right-hand  side  of  (1), 
one  may  note  that  the  N  principal  roots  of  z**-e^b^  are 


i 

e 


bT-2ir 

N 


i- 


bT+27r 

N 


j 

e 


bT-Air 

N 


1  T/N 

Hence  the  partial  fraction  expansion  of  — C  (z)  may  be  written  as 

z 


and 


±CT/V) 

z 


n2  A  +  jB 

l  -  "  ° 


n=nl  2- 


jw  T/N 


Z  -  E 


+  [terms  due  to  modes  of  (GM)^N] 


(6) 


where 


,  ,  2irn 
u>  ”  b  +  — — 
n  T 


n^  *  — ^(N-l) 


n2  -  y(N-l) 


if  N  *  odd 


2 


I- 


if  N*  even 


A»  +  - 


tan 1/11  z”-1 


z*e 


jut  T/N 
n 


z*e 


junT/N 


^(GM(s))T/N 


s-ju 


(7) 


The  steady-state  response  to  e^bt  can  be  written  by  inspection  of  (6)  as 


Ic88<t>] 


T/N 


T  r2  ju>  t*lT/N 


(8) 
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Since  the  system  of  Fig.  1  is  linear  the  response  to  the  imaginary  part 
of  is  the  imaginary  part  of  the  response  to  e^C.  Hence  the 

steady-state  response  to  slnbt  is 


Icsa(t)] 


T/N 


Tp 

I  >  (A  sin  u  t  +  B  cosw  t) 
I  L  n  n  n  n 

ln*ni 


T/N 


(9) 


where  and  Br  are  determined  by  (7).  Note  the  simplicity  of  the 
derivation  as  compared  to  that  of  Reference  [1].  The  difference  between 
(9)  and  the  corresponding  expression  (Eq.  (11)  of  Reference  [1])  is 
that,  in  the  right-hand  side  of  (9),  the  N  spectral  components  are 
selected  alternately  between  negative  (subalias)  and  positive  aliases, 
while  only  positive  aliases  are  used  in  Eq.  (11)  of  Reference  [1],  It 
is  worth  noting  that  while  both  results  are  correct  for  the  representa¬ 
tion  of  sampled  steady-state  output,  difficulty  arises  (and  was  over¬ 
looked)  in  considering  the  limit  case  of  N-*-»  in  Eq.  (11)  of  Reference 
[1].  By  letting  N -*■<*>  in  (9),  this  difficulty  is  circumvented.  The 
following  development  demonstrates  this  point.  The  correct  expression 
resulting  from  taking  N-»-«  in  (9)  is  presented  first. 

Letting  N-*-®  in  (9)  gives 

00 

c  (t)  ■  Y  (A  sin  w  t  +  B  cos  w  t)  (10) 

ss  L  n  n  n  n 

n™— ® 

where,  on  account  of  (7), 

An  +  Bn  ■  711" 


b  lim  ! 

T  N-**»  k- 


GM(s+j~^) 


s-jw 


(ID 
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N  N 

As  long  as  -y£n^  in  in^  and  N-*-*,  only  the  k«0  term  contributes  to 
the  infinite  sum  in  the  right-hand  side  of  (11)  since  GM(s)  is  always 
a  low-pass  filter.  Hence,  as  N-*», 

An  +  jBn  -  IfflO-n)  U2> 

However,  the  above  limit  case  is  not  true  if  n*N-i,  for  finite  integer 
i.  In  that  case, 

A  +  jB  -  ^  lim  l  GM(s  -  j~~~) 
n  n  i  ,,  ,  x 

N-**3  k"-°° 

-  4  lim  l  GM[j(b  + 

N-*°  k=— “ 

-  i  lim  l  GM[j(b  -  J2-(-kf--)N]  (13) 

N-x”  k=>-“ 

Now  as  N  -*■“>,  only  the  k-1  term  contributes  to  the  infinite  sum  in  the 
right-hand  side  of  (13)  since  GM(s)  is  always  a  low-pass  filter.  Thus 
for  n*N-i,  N -*■«>,  finite  i, 

An  +  jBn  -  YGMlj(b  "  (14) 

which  is  a  subalias  component,  is  not  negligible,  and  is  not  represented 

by  (12).  Since  the  summation  in  Eq.  (11)  of  Reference  [1]  runs  between 

n*0  and  n“N-l,  letting  results  in  an  A  +1B  which  cannot  be 

n  n 

represented  by  (12),  which  is  Eq.  (18)  of  Reference  [1). 

A  numerical  example  demonstrating  the  representation  of  the 


spectral  components  of  the  continuous  output  and  the  sampled  output  on 
a  Bode-plot  and  an  even-scaled  frequency  response  plot,  respectively. 


has  been  given  in  Reference  [1],  The  modification  and  correction  made 
in  this  paper  requires  the  inclusion  of  subaliases  or  negative  aliases 
on  the  frequency  response  plots  in  accordance  with  (8)  and  (10).  To 
avoid  duplication,  numerical  examples  will  not  be  given  here.  The 
reader  is  urged  to  apply  this  modification  to  the  example  given  in 
Reference  [1]. 

III.  DIRECT  DERIVATION  OF  THE  FREQUENCY  SPECTRUM  OF  THE  CONTINUOUS 
OUTPUT 


In  the  above  development,  the  continuous  output  is  treated  as  a 
limit  case  of  sampled  output  with  an  arbitrarily  large  sampling  fre¬ 
quency,  and  the  infinite  frequency  spectrum  of  the  continuous  output  is 
given  by  (10)  with  the  spectral  components  An+jBn  given  by  (12) .  In 
this  section,  it  will  be  shown  that  the  frequency  spectrum  of  the  con¬ 
tinuous  output  of  a  discretely  excited  continuous  system  can  also  be 
derived  directly  without  invoking  the  sampled  output.  Thus  it  confirms 
the  result  of  the  previous  section. 

Again  consider  the  system  of  Fig.  1,  where 

C(s)  -  G(s)M(s)RT(s)  (15) 

Let 

gm(t)  -  L"1[G(s)M(s)]  (16) 

Let  the  input  be  e^C,  with  t  extending  from  -•  to  «.  Then  the  steady- 
state  response  of  the  system  may  be  written  as 

cas(t)  -  [eJbt]T*ihl(c>  (17) 
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where  the  *  denotes  convolution,  and 


2irn 

cjbt  i  S(t-nT)  -  itJbt  j  S 
n=-»  n=—» 


(18) 


The  last  identity  can  be  found  in  textbooks  on  Fourier  series  and 
Fourier  transform  [2].  Fourier  transformation  of  (18)  gives 


f[(eJ1,C}T]  -  I 


T  L  6  (w  -  b  -  ~) 
n=-°° 


(19) 


Fourier  transform  of  (17)  is 


CsS(j“)  =  T  ^  6(w'b 

n=-® 


2ml 


)G(  jo))M(  jw) 


(20) 


since  g  (t)  is  a  causal  function.  Now  the  steady-state  output  is 


m 


obtained  by  inverse  transforming  (20) : 


1  r  2irn  + 

css(t)  "  T  I  GMlKb  +^)]  e  T 

n=-» 


(21) 


ry 

Since  a  b  +  — substituting  (12)  into  (21)  gives 


®  ju  t 

c  (t)  -  y  (A  +  jB  )  e  n 
ss  L  n  J  n 

n"  - 


(22) 


as  the  steady-state  response  to  e^bt.  Note  that  (22)  is  the  limit  case 
of  (8)  when  N-+-*.  Again,  by  virtue  of  the  linearity  of  the  system, 

(22)  implies  that  the  steady-state  response  to  sin  bt  is  given  by  (10), 
with  An+jBn  given  by  (12).  In  summary,  the  analogy  between  the  fre¬ 
quency  response  of  a  continuous  system  and  that  of  a  discretely  excited 
continuous  system  is  readily  demonstrated  by  comparing  Fig.  2  and  Fig.  3. 
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Note  especially  that  if  all  the  variables  are  represented  in  the 
s-domain,  a  unified  concept  of  frequency  response  may  be  readily  applied 
to  both  continuously  excited  and  discretely  excited  continuous  systems. 


i 


IV.  CLOSED-LOOP  DIGITAL  CONTROL  SYSTEMS 


It  is  readily  seen  that  if,  in  Fig.  2, 


F(s) 


Gm[i  +  G^(G2GM)T 


1-1 


(23) 


then  Fig.  2  represents  the  closed-loop  system  of  Fig.  4.  Since 
wn=b  +  -^p-  and  the  Laplace  transform  of  a  sampled  function  is  periodic 
in  s  with  period  j^p,  substituting  s=*jwn  into  F(s)  gives 


F(>n) 


GM(ja»n)|[: 


1  +  G*(G2GM) 


•(24) 


The  concept  of  frequency  response  set  forth  in  Fig.  2  implies  that,  for 
the  continuous  output  of  the  closed-loop  system  of  Fig.  4, 

A„  +  ‘  +  GI<G2GM>T]  ^  <«) 

v  Js*jb 


For  the  sampled  output  with  sampling  period  T/N,  the  spectral  components 
are  given  by  substituting  the  right-hand  side  of  (23)  for  GM  in  (7). 

Thus 


(26) 


Applying  multirate  sampling  theory  (Appendix  of  [1])  and  property  of 
periodicity  to  (26)  gives 
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A»  +  i».  ■  [s<OT)T/N]3.j(Jii  (t1  +  *V*f*}^  <”> 

Thus  the  spectral  representation  of  a  closed-loop  discretely  excited 
continuous  system  is  readily  obtained  via  the  unified  concept  of  fre¬ 
quency  response. 


V.  CONCLUSIONS 

This  paper  corrects  the  mistake  made  in  a  previous  paper  and 

further  advances  the  concept  of  frequency  response  of  discretely  excited 

continuous  systems.  If  a  continuous  system  has  transfer  function  F(s) 

and  its  input  is  R(s),  then  its  steady-state  response  to  r(t)=e^bt  is 

c  (t)=F( jb)e^bt  (Fig.  3).  If  R(s)  is  sampled  and  then  input  into  F(s), 
ss 

then  the  input  to  the  continuous  system  F(s)  is  (Fig.  2) 

RT(s)  «  Y  l  R<s+j^>  (28) 

n=-<*>  1 

and  the  steady-state  response  of  the  discretely  excited  system  to 
r(t)*eJbt  is 

css(t)  “  T  l  Ftj(b  +  e  (29) 

n“-»' 

The  analogy  between  the  frequency  response  of  a  continuous  system  and 
that  of  a  discretely  excited  continuous  system  is  interesting  in  the 
sense  that  it  enables  one  to  write  the  frequency  response  of  the  latter 
on  the  basis  of  what  is  already  known  about  the  former.  This  unified 
concept  of  frequency  response  is  independent  of  system  conf igurations, 
as  F(s)  may  represent  closed-loop  as  well  as  open- loop  systems. 
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n  =  -co 


Fig.  2  Frequency  Response  of  Discretely  Excited  Continuous 
System 


If  r  (t)  =  «ibf  Then  Css  (»)  *  F  (jb)  «ibl 


Fig.  3  Frequency  Response  of  Continuous  System 
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p  Digital  Control  System 


REFERENCES  FOR  APPENDIX  D 


[1]  Whitbeck,  Didaleusky  and  Hofmann:  "Frequency  Response  of  Digitally 
Controlled  Systems,"  J.  of  Guidance  and  Control,  Vol.  4,  No.  4, 

July -August  1981,  pp.  423-427. 

[2]  McGillem,  C.  D.  and  Cooper,  G.  R.:  Continuous  and  Discrete  Signal 
Analysis,  Holt,  Rinehart  and  Winston,  1974,  p.  158. 

[3]  Ragazzini,  J.  R.  and  Franklin,  G.  F.:  Sampled-data  Control  Systems, 
McGraw-Hill,  1958. 


80 


END 

DATE 

FILMED 


